Independence and interdependence in the nest-site choice by honeybee swarms: 
agent-based models, analytical approaches and pattern formation 
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In a recent paper List, Elsholtz and Seeley [List et al. 20091 ] have devised an agent-based model of 
tlie the nest-choice dynamics in swarms of honeybees, and have concluded that both interdependence 
and independence are needed for the bees to reach a consensus on the best nest site. We here present 
a simplified version of the model which can be treated analytically with the tools of statistical physics 
and which largely has the same features as the original dynamics. Based on our analytical approaches 
it is possible to characterize the co-ordination outcome exactly on the deterministic level, and to 
a good approximation if stochastic effects are taken into account, reducing the need for computer 
simulations on the agent-based level. In the second part of the paper we present a spatial extension, 
and show that transient non-trivial patterns emerge, before consensus is reached. Approaches in 
terms of Langevin equations for continuous field variables are discussed. 
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I. INTRODUCTION 



In a recent paper List et al. 2009l | List, Elsholtz 



and Seeley (LES for short from now on) presented an 
agent-based model of the decision making processes 
in the nest-site choice of honeybee swarms, and their 
paper has attracted a significant amount of atten- 
tion in the public media (see e.g. The Economist of 
13 February 2009). LES focus on the interplay be- 
tween the interdependence among bees and their in- 
dependence in assessing the quality of potential nest 
sites, and conclude that, assuming a degree of inde- 
pendence, bees will generally choose the best nest site 
for a wide range of non-extremal model parameters, 
and secondly that both independence and interdepen- 
dence are necessary for the swarm to reach a reliable 
consensus on the best possible nest site. The precise 
notions of 'independence' and 'interdependence' will 
become clear below, but in essence LES find that, in 
order to identify the best nest site, the bees in their 
model need to make independent assessments of the 
true quality of nest sites instead of blindly following 
the choices advertised by fellow bees. At the same 
time interdependence, i.e. direct communication be- 
tween bees, is required as well in order to promote a 
quick and accurate convergence towards a consensus 
on the best choice. 

While the work of LES is mainly based on com- 
puter simulations, we here discuss a modification of 
their agent-based dynamics, which can be addressed 
using the tools and techniques of statistical mechan- 
ics. Such analytical solutions mean that computer 
simulations on the agent-based level could in princi- 
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pie be dispensed with entirely, although admittedly 
our resulting equations still require numerical solu- 
tion. Still, if there are k nest sites in the model, then 
on the deterministic level, the solution of the reduced 
model comes out as a coupled set of k quadratic equa- 
tions for the number of bees advertising the different 
sites, which is a significant reduction in complexity of 
the problem. LES choose fc = 5 in their analysis, so 
that we are able to reduce the deterministic features 
of their model to a very small set of equations. Of 
course, stochasticity is important in the model, but 
given the deterministic solution, fluctuations about 
the deterministic limit can be characterized analyti- 
cally as well to a good approximation, as we will de- 
scribe below. The main aim of this current paper is 
to present the modified version of LES's model, and 
to show how it can be addressed with the tools of 
statistical physics. In the second part of o ur work we 
then extend the model of [List et al. 2009| to include 
space, and study the pattern formation dynamics of 
this extended model. 

The natural process modeled by LES is the follow- 
ing: once a year once a colony of bees has reached 
a certain size, a 'search committee' consisting of sev- 
eral hundred bees swarm out to identify suitable nest 
sites, and once an agreement has been reached the 
original colony divides, and the queen leaves with 
about two third s of the bees for the new nest (see 
[List et al. 2009j | and references therein). The pro- 
cess LES are interested in is the actual decision mak- 
ing. The bees who fiy out to inspect potential new 
sites return to the original colony and perform dances 
to advertise the new sites they have inspected. The 
duration of the dance here is a measure for the qual- 
ity of a newly inspected site, the better the site, the 
longer the dance. Fellow bees who observe these 
dances might then themselves start dancing for that 
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particular site as well, and this may happen with or 
without an independent inspection of the actual site, 
i.e. with or without an independent assessment of 
the site's quality. This may then lead to other bees 
initiating dances (again with or without independent 
assessment of the site's quality), and finally, if all goes 
well, a general consensus on one of the nest sites will 
be reached. The question LES address in their model 
is whether this consensus will be on the best possible 
site, or whether a sub-optimal choice can result. 

In order to model communication and interdepen- 
dence between bees LES assume that the rate with 
which a given bee starts dancing for a particular site 
depends on an a-priori rate of independently discov- 
ering that site, and, weighted by an 'interdependence 
coefficient', the fraction of other bees dancing for this 
site. If the interdependence between bees is strong 
then the latter factor carries a large weight, and the 
rate with which bees start dancing for a given site 
is mostly determined by the fraction of bees already 
dancing for that site. If interdependence is weak 
dances for all sites are commenced essentially with 
equal rates, independently of how many bees are al- 
ready advertising the different sites. The second com- 
ponent of the model is the time a dance lasts once it 
has been begun. In the model by LES, the dura- 
tion of the dance is determined by a combination of 
the actual (or perceived) quality of the site, and a 
uniform site-independent contribution. In the case of 
strong independence, the bees are assumed to inspect 
the site before dancing for it, and the duration of the 
dance will be proportional to the perceived quality 
of the site. In the case of weak independent quality 
assessment the duration of the dance for a given site 
does not depend much on the site's quality. 

The key question LES address in their study is 
whether or not the swarm of bees will eventually 
reach a consensus on what the best nest-site might be, 
and whether this consensus choice then is the actual 
best site. The outcome of the decision dynamics here 
depends on the independence and interdependence 
parameters introduced above, and, to a lesser extent, 
on the accuracy with which bees assess the quality of 
a site if they actually inspect it. I.e. a bee may decide 
to visit a site before dancing for it, but the perceived 
quality may be a random function of the actual qual- 
ity, so that noise may blur the bee's assessment of 
the site. In this paper we complement the simulation 
study of LES by analytical computations, allowing for 
a more extensive characterization of the model. We 
also present extensions towards a spatial model, and 
discuss how patterns may emerge at transient times 
if the model is considered on a two-dimensional regu- 
lar lattice. The paper is organized as follows: In Sec. 
im we will first re- iterate the definitions of the model 
by LES, and then introduce the simplified dynamics, 
which we will then address analytically within a mas- 
ter equation approach in Sec. IIIII We here formulate 



deterministic mean-field equations and also address 
first-order finite-size corrections. In Sec. IIVI we then 
introduce a spatial version of the reduced model, and 
show in computer simulations that non-trivial tran- 
sient patterns may emerge if model parameters are 
chosen suitably. The study of the spatial model is 
complemented by the numerical integration of a set 
of partial differential equations for a coarse-grained 
dynamics. In Sec. |V] we finally summarize our work 
and point to future directions which may be of inter- 
est. 



II. MODEL DEFINITIONS 

A. Original model by List et al 

In the original model by LES there are N bees, 
labelled i = 1,2, . . . , N, and k potential nest sites, 
we label them by a = 1, . . . ,k. Each nest site has 
an intrinsic quality > 0. These 'actual' quali- 
ties of the different nest sites are fixed from the be- 
ginning and do not change over time. We always 
choose qi < q2 < ■ ■ ■ < qk, i.e. site number k 
is the best site, and site number 1 the worst. Ad- 
ditionally, for each potential nest site a, there is 
an a-priori probability Tr^ of finding this particu- 
lar site. These rates of finding the different nests 
will be taken to fulfill X]a=i ^ 1' ^"^"^ ^'^^^ 
set ttq = 1 — > .„_] TTn,. T hroughout this paper we 



will follow [List et al. 2009l | and assume that the it a, 
a = 1, . . . ,k, are identical for all nests, but an exten- 
sion of the model to the more general case is straight- 
forward. 

The decision and dancing dynamics within the 
model then proceeds in discrete time steps, t = 
1,2,3,..., and at each time step each bee can ei- 
ther be dancing for a particular site a € {1, . . . ,k}, 
or not be dancing. The not-dancing state here de- 
scribes bees that may have flown out to inspect a 
site, but have not yet started dancing for it, as well 
as bees wh o are observing other bees, or bees that 



are resting List et al. 2009l |. We will sometimes re- 



fer to these bees as 'passive'. Following LES we write 
Sj(i) = a e {1, . . . , fc} if bee number i is dancing for 
site a at time t, and Si{t) = 0, if she is not dancing at 
time t. In LES's model, the state of a bee is character- 
ized by an additional variable, di{t) > 0, indicating 
the remaining duration of bee I's dance at time t. I.e., 
up to rounding to a near integer, if di{t) > 1, then 
bee i will be dancing for another di{t) time steps. In 
the non-dancing state (i.e. for a bee i with Si{t) = 0) 
di{t) has no relevance. More specifically, the dynam- 
ics proceeds according to the following algorithm: 

1. At time t ~ initialize the {si{t),di{t)}, for 
example set Si{t = 0) = 0, and di{t) = for all 
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2. At time t compute the fraction of passive bees 
and of those dancing for the different sites {a} 

at this time, fa{t) = N-^J2f^iSa,s,{t) S [0,1] 
(a = 0,...,k). S here stands for the Kro- 
necker delta, i.e. 6a[3 = 1 if a = /3, and 
Saf3 = otherwise. Then compute pa(t + 1) = 
(1 - A)7r„ + Xfait) {a = 0,...,k). Note that 
^a=oPo:i^ + 1) = 1 by construction, due to the 
normahzation of the tt^ and of the fa- The 
variable A G [0, 1] is a model parameter, and its 
meaning will be explained below. 

3. In order to update the states of the bees dis- 
tinguish between bees that are not dancing at 
time t and those which are dancing at time t. 
Update all bees in parallel dynamics, i.e. for all 
i = 1, . . . , N perform the following: 

3a) If bee i is dancing at time t, then check 
whether di{t) > 1. If this is the case then 
the bee keeps dancing for the site she is 
currently dancing for, i.e. set Si{t + 1) = 
Si (t) and reduce the remaining duration of 
the dance by one unit: di{t+l) = di{t) — l. 
If the bee is dancing, but di{t) < 1, then 
set Si{t+1) = and di{t + l) = (the bee 
stops dancing). 

3b) If, at time t, bee i is not dancing {si{t) = 
0), then with probability pa{t + 1) set 
Si{t + 1) = a {a = l,...,fc), i.e. with 
probability pi(i + l) the bee starts dancing 
for site a — 1, with probability p2(t -\- 1) 
she starts dancing for site a — 2 and so 
on. With probability pQ{t + 1) she will re- 
main in the non-dancing state. If a dance 
is commenced then it remains to specify 
the duration of the dance bee i has just 
started for site a. Here draw a random 
number T from a standard Gaussian dis- 
tribution (mean zero, unit variance). Then 
set 

{Qa exp((TT) with prob. 1 — /i 
K explaT) with prob. /i 

4. Iterate, i.e. go to 2. 

LES here set K ~ qk, i.e. K is the quality of the best 
site. The variables A € [0, 1] and 1 — /i S [0, 1] are the 
interdependence and independence parameters men- 
tioned in the introduction. If A is close to one, then 
the interdependence between bees is strong, i.e. the 
rate of commencing a dance for a given site is mostly 
determined by the number of bees already dancing 
for this site. At small A the dancing activity of other 
bees is mostly irrelevant for a given individual to take 
up a dance. The independence in the assessment of 
the nest-site qualities is parametrized by fi, more ac- 
curately by 1 — /i in the notation of LES which we 




FIG. 1: Original model by List et al: Fraction fa of 
bees dancing for each of the sites, a = 1,2..., 5, in 
the stationary state, along with the fraction of bees who 
are not dancing, /o, versus A at fixed /i = 0. Mark- 
ers from simulations (TV — 200, run for 1000 steps, av- 
eraged over 10 runs). Model parameters are fc = 5, 
(go, <7i, 92, (73,94,(75) = (0.1,3,5,7,9,10), a = 0.2, and 
ttq — 0.05 for all a = 1, . . . , 5. 



will follow here. For large 1 — ^J- the duration of a 
dance for site a strongly depends on the quality Qa of 
that site (independent assessment), whereas for small 
1 — /i, the actual quality of the different sites is mostly 
irrelevant for the duration of the dances. The vari- 
able (T > is a further model parameter and rep- 
resents the reliability of the bees' assessment of the 
qualities of the different sites. As seen in the above 
update rules the perceived quality of site a is given 
by Qa exp((7T), where T is a standard Gaussian ran- 
dom variable. Thus if cr = 0, the perceived quality 
is always identical to the actual quality Qa- Increas- 
ing a introduces more and more uncertainty into this 
perception. 

The dynamics as defined above is easily imple- 
mented in a computer simulation, and LES have pre- 
sented an extensive numerical analysis of the model. 
For completeness we present some simulation results 
in Fig. [Tl where we show the fraction of bees danc- 
ing for the different sites in a situation where there 
are five potential nest-sites, and with strong indepen- 
dence (1 — /i = 1), and varying interdependence. As 
seen in the figure, the best nest-site tends to be the 
most populated one on average, but convergence to 
the best nest-site is strong only at moderate to high 
interdependence, but not so much at low values of A. 
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B. Simplified model 

If one introduces state variables Xi{t) = 
{si{t),di{t)), as LES do, then the above dynamics 
defines a Markov chain in the state space spanned by 
the {xi}, i = 1,...,N, and a master equation ap- 
proach is feasible in principle. However, as we will 
discuss in the next section, the model can actually be 
reduced first, without affecting its key features and 
behavior much, to yield a dynamics which is Marko- 
vian in the space of the {si{t)} alone. This simplifies 
the analysis considerably, and as we will show next, 
the 'memory' variables di{t) can be dispensed with 
entirely. 

The simplified model operates purely in the space 
of the {si{t)}, i.e. the state of the system at 
any given time is fully determined by the variables 
si(<), . . . , SAr(t). In essence the simplification con- 
sists in the following step: if in the original model 
a bee begins to dance for a given site at a given time 
t, the duration di{t) of this dance is determined at 
the beginning of the dance, and then, up to round- 
ing to a near integer, at time step t + di(t) the bee 
stops dancing for this given site. For /i < 1, the du- 
ration of the dance depends on the quality of the site 
danced for, and is also subject to some uncertainty 
as explained above. But crucially, the duration of the 
dance is set at the beginning of the dance, the end 
of the dance then occurs at a deterministic moment 
in time, once the duration has been chosen. The du- 
ration of the dance may have a stochastic element 
to it as explained above, but its value is fixed and 
known when the dance starts. In the reduced model 
we replace this with a random process, i.e. a bee may 
start dancing for a given site a at time t, and then 
at each subsequent time step the bee ceases dancing 
with a certain probability and continues dancing 
with probability 1 — Tq, . The rate Tq, with which the 
dance is stopped at each step may here depend on 
the quality of the site, for example a high-quality site 
a will have a lower value of Ta than a site /3 of a 
lesser quality. The probability that a dance for a 
given site a lasts for precisely £ time-steps after the 
dance has begun is then given by {l — raY^-^ra, i.e. it 
follows a geometric distribution of mean l/ra- This 
may be slightly unrealistic in the sense that the prob- 
ability of stopping after exactly £ rounds is a mono- 
tonically decreasing function of £, whereas the distri- 
bution of dance durations used by LES has a well- 
defined maximum, but this drawback is, in our view, 
largely outweighed by the analytical simplifications 
it brings with it. As we will show below this sim- 
plification does not affect the qualitative behavior of 
the model much. A further restriction of the reduced 
model is the fact that there is no clear analogue of 
the uncertainty parameter a in the present setup, i.e. 
once Tq, has been chosen, both the mean and variance 
of the distribution of durations of dances for site a 



are fixed. We will discuss possibilities to include such 
uncertainty towards the end of this paper. 

The precise algorithm of the simplified model is the 
following, the definitions and conventions regarding 
the {tTq} and the model parameters A and fi are the 
same as in the original model: 

1. At time t = initialize the {si{t — 0)} for all i, 
for example set all bees in a non-dancing state, 
or initialize the Si{t) at random. 

2. At time t compute the fraction fa{t) € [0, 1] of 
bees dancing for site a at time t (a ~ 0, . . . , k), 
and compute Pa{t + 1) = (1 — A)7rQ + Xfa{t). 
Note that < ^^(i + 1) < 1 for all a, and 
Z]a=oP"(* + 1) = 1 by definition. 

3. In order to update the states of the bees dis- 
tinguish between bees that are not dancing at 
time t and those which are dancing at time t. 
Update all bees in parallel dynamics, i.e. for all 
i — 1, . . . , N perform the following: 

3a) If bee i is dancing for site a at time t (i.e. 
if Si{t) = a > 0), then with probability Tq, 
set Si -I- 1) = (bee stops dancing), and 
with probability 1 — set Si{t + I) = a 
(bee keeps dancing for site a). 

3b) If, at time bee i is not dancing {si{t) = 
0), then with probability Pa{t + 1) set 
Si{t + 1) — a, a — 0, 1, . . . , /c, i.e. with 
probability pi{t + 1) the bee starts danc- 
ing for site a = 1, with probability p2(^+l) 
she starts dancing for site a = 2 and so on. 
With probability po{t + l) she will remain 
in the passive state a = 0. 

4. Iterate, i.e. go to 2. 

It remains to specify the rates with which dances are 
terminated, i.e. to define the {tq,}. We here choose 



K 



a 



List et al. 2009l |. 



where K = as in 

the rate Tq, is a decreasing function of 



,...,A:, (1) 

For any /i < 1 
i.e. the 



higher a site's quality, the lower the rate at which 
dances for it are stopped. The pre-factor go (which 
needs to be non-zero) here ensures that < < 1 
and also that the {tq,} do not depend on the overall 
scale of the nest-site qualities. Equivalently, go, to- 
gether with the rate at which dances are commenced, 
defines the time-scale of the model. A different scal- 
ing of time could be chosen by using 759 as a pre- 
factor in ll]), with < 7 < 1. We note that for /i = 
(maximal independence) one has — qo/qa, i-c. 
is (up to the pre-factor) given by the inverse quality 
of the nest-site a, so that the mean duration of dances 
is proportional to the nest-site quality. For fJ. — 1 (no 
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FIG. 2: Temporal decision dynamics of the group of bees 
at = 0, A = 0.8. We show the fraction fa{t) of bees 
dancing for each of the sites, a = 1,2,. ..,5, along with 
the fraction of bees who are not dancing. Solid lines are 
from the deterministic dynamics, Eqs. (0, markers from 
simulations started from an all-passive state (A'' = 400 
bees, averaged over 250 runs). 



FIG. 3: Fraction /* of bees dancing for each of the sites, 
Q = 1, 2 . . . , 5, in the long-run, along with the fraction of 
bees who are not dancing, versus A at fixed fi = 0. Solid 
lines are from the deterministic dynamics, markers from 
simulations (TV = 400, run for 20000 steps, averaged over 
50 runs). 



independent assessment of the nest-site quality), one 
has = qo/K irrespective of a G {1, . . . , fc}, so that 
the mean duration of dances for a given site does not 
depend on its quahty at all. It is here worth point- 
ing out that if the initial state is chosen to be the 
all-passive state, then for A = 1 no bee ever starts 
dancing. An analogous effect is present in the orig- 
inal model by LES. We will therefore generally ex- 
clude the case of A = 1 from the analysis, and choose 
A large, but strictly smaller than unity, if we want to 
model strong interdependence. Alternatively one can 
choose random initial conditions. 



III. ANALYTICAL SOLUTION AND TEST 
AGAINST SIMULATIONS 



Master equation and dynamics in the 
deterministic limit 



In order to analyze the stochastic process de- 
fined by the dynamics of the model, we will write 
na{t) for the number of bees dancing for each of 
the sites at time t. We here explicitly include the 
state a = 0, i.e. bees not dancing at time t. One 

then has X)q=o"'"(0 — for all times t. For 
later convenience we introduce the vector n{t) — 
{no{t),ni{t), . . . ,nk{t)). The state of the system at 
time t is therefore fully determined by n(t). The re- 
duced model defines a Markov process for n(i), de- 
scribed by the following master equ ation (see text- 
books such as [van Kampen 199^ or Gardiner 2009l | 



for details on the master equation formalism) 



^Pn = ^P„_e„+eoT+(n-e„-t-eo)-P„^T+(n) 



^n+e„-eoT^a (n + e„ - Gq) - P„ ^ (n) 



(2) 



for the probability Pn{t) of finding the system in state 
n at time t. T+(n) here stands for the probabil- 
ity that a randomly chosen bee is inactive and starts 
dancing for site a £ {I, . . . ,k} at time t, given that 
the system is in state n. Similarly, ^^(n) is the prob- 
ability for a randomly chosen individual to be dancing 
for site a and for it to cease dancing in the subse- 
quent time step. We write e^, a = 0,1,..., A; for 
the k + 1 unit vectors, i.e. one has (e^)^ — 5ap for 
a, P G {0, 1, . . . ,k}. The transition rates are given by 



and 



r-(n) 



no 
N 



N 



90 



ii 

K 



(3) 



(4) 



respectively, where a — l,...,fc. Technically 
speaking, this master equation is more appropriate 
for so-called random sequential updating, in which 
the states of the N bees are updated one-by-one 
(and each update would then take a time-step of 
1/iV, so that each bee is updated on average once 
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per unit time) , see [van Kampen 1993 . iRisken 19961 . 
iGardiner 2009l | for details. The pre- factors Ua/N in 
the rates here in fact represent the probabilities 
that a randomly chosen bee is in state a. We have 
verified in numerical simulations that there are no sig- 
nificant differences between the stationary states of 
the model with parallel and sequential updating re- 
spectively, and as we will see below theoretical predic- 
tions from our master equation approach agree well 
with simulations using parallel update rules. 

Introducing continuous variables fa = ria/N for 
a = 0, 1, . . . , fc it is then straightforward to derive a 
set of deterministic ordinary differential equations in 
the limit of an infinite system size, N —f oo. This 
may for example be done by considering 



solution of the fc coupled quadratic equations. 



(n(i))=^nP„(i) 



(5) 



and noting that 

k 

-t-^[(eo-e„)T-((n))], (6) 

Q = l 

where we have used a deterministic approximation to 
write (r±(n)) =T±((n)). This gives 

Ut) = (I-pW) [(1 - A)^„ + A/„(t)]-r„/„(t), (7) 

for a = 1, . . . , fc, and where p{t) = ^2^=1 fai^)- The 
fraction of bees not dancing at time t is given by 

B. Comparison with simulations 

Equations ([7]) are easily integrated numerically for 
any fixed choice of the model parameters, and we 
show a specific example in Fig. [21 We her e choose 
the model parameters as in [List et al. 2009l | , specifi- 
cally we use fc = 5, (go, ft, ■•■,<?5) = (0.1,3,5,7,9, 10), 
K = (75, and tti = ... = 775 = 0.05, as well as 
ttq = 0.75. Our definitions of the transition rates re- 
quire (Jo > and so we choose a small, but non-zero 
value for go- As mentioned above go sets the overall 
time-scale of the model, and is essentially arbitrary, 
so that the times indicated in Fig. [5] do not necessar- 
ily have a direct biological interpretation in terms of 
e.g. hours, days or the total number of nest inspec- 
tion flights done by the population. As seen in the 
figure the theoretical description of Eqs. ([7]) agrees 
near perfectly with results from simulations, and re- 
produces even non-monotonic trajectories faithfully. 
The stationary state of the system may therefore be 
determined as a fixed point of Eqs. (O, i.e. as the 




A)7r„ + Xa (8) 



for a = 1, . . . , fc, and where we have introduced the 
notation /* to denote the fixed-points of ([7]). These 
fixed-point equations may in principle be solved nu- 
merically, or alternatively, one integrates Eqs. ([7]) to 
asymptotic times. We have here followed the latter 
approach. It is interesting to note that Eqs. ^ can 
be re-arranged to give 



1 - p* 



l-A vr, 



,fc. (9) 



Given that all sites have the same a-priori probabil- 
ity of being found, i.e. that tTq does not depend on 
a e {1, . . . , fc}, and that secondly the rates of dances 
being terminated decrease with an increasing qual- 
ity of the sites, Va < rp for a > /3, Eqs. (O re- 
veal immediately, that a site is the more populated 
on the deterministic level, the higher its quality, i.e. 
A* < /I < • • • < /fc, in-line with the results of LES. 
We also note that for A = 1 the only possible fixed 
point is one at which = for all a = 1, . . . , fc. 

To test our theoretical predictions regarding the 
stationary state we show the asymptotic values of 
/o,/i,...,/5 as a function of the model parameter 
A in Fig. [3l and again simulations confirm that our 
mean-field description accurately reproduces the be- 
havior of the agent-based model in a wide range of 
model parameters. Figure [3] also reveals that the be- 
havior of the simplified model is very similar to that 
of the origin al model by LES (see Hypothesis 1 and 
Figs. 1-3 of [List et al. 2m^ \. if the bees assess the 
quality of the sites independently before dancing for a 
given site, i.e. if ^ = 0, then the best site is retrieved 
successfully over a wide range of val ues for the pararn - 
eter A. In-line with the findings of List et al. 2009l |. 



an increased amount of interdependence between the 
bees (high values of A) increases the chances of find- 
ing the best nest site, as indicated by the proportion 
of bees, /g , dancing for the best site a — b. 

Having established the correctness of our theoreti- 
cal approach, we can now push the analysis further, 
and assess the decision-dynamics of the bees through- 
out the parameter space defined by A and p,. To 
this end we have integrated Eqs. (O to large times 
[t = 10^) for approximately 1.6 • 10^ different combi- 
nations of A and /i, equally spaced in the unit square, 
and have computed the difference A = /g — /| be- 
tween the fraction of bees dancing for the best site 
and the fraction of bees dancing for the second-best 
site. A therefore ranges in the interval [0,1], and 
high values of A indicate that the best site [a = 5) 
is identified with a high accuracy by the population 
of bees. Results are shown in Fig. [H and confirm 
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FIG. 4: Quality of the retrieval of the best site, as mea- 
sured by A (see main text), throughout the parameter 
space. High values of A indicate a high accuracy in iden- 
tifying the best nest site. Results are from an integration 
of Eqs. (0, using a first-order Euler-forward scheme with 
a time-stepping of = 0.1, integrated up to time t — 10^. 



that a reliable convergence to the best site can only 
be achieved in the presence of both, a sufficient de- 
gree of independence in the assessment of nest-sites 
and interdependence between the different individ- 
uals. This is a confirmatio n of the second main hy- 
pothesis of List et al. 2009l | . The best retrieval of the 
optimal site is found at low values of jji (high indepen- 
dence) and large A (high interdependence), fn order 
to check for convergence we have also analyzed data 
from an integration of Eqs. ([7]) only up to t = 10'*. 
The location of the different lines in Fig. [D remains 
unchanged at this lower time, so that we conclude 
that the stationary state has already been reached 
much earlier than at t = 10^. 



C. Size-limitations of best-site retrieval 

The deterministic dynamics defined by Eqs. ([7]) ap- 
proaches a stationary state in which the proportion 
of bees populating the best site is higher than that of 
any other site. The description in terms of these de- 
terministic equations is however only valid in the for- 
mal limit of an infinite population of bees, i.e. if iV ^ 
oo. In small populations stochastic effects may in- 
stead dominate the decision making, and as observed 
by LES convergence to a sub-optimal state may oc- 
cur. To characterize these finite-size effects one may 
look at a system at a fixed size N and at fixed model 
parameters A and ^l. At any given time, the state 
of the system is described by the (re-scaled) pop- 



ulation vector {ni{t)/N,n2{t)/N, . . . ,nk{t)/N), and 
the dominating choice at any given time may well 
be different from the optimal site, i.e. it may well 
be that the most-populate site at time t is not ac- 
tually the best site a = k. This is the case when- 
ever there is an index a G {l,...,fc — 1} such that 
na{t) > nk{t). It is therefore instructive to look 
at the joint and marginal distributions of the values 
{ni{t)/N,n2{t)/N, . . . ,nk{t)/N) sampled over time 
and from different simulation runs. We first ad- 
dress the marginals, and the resulting histograms are 
shown for a particular choice for A and /x and for a 
relatively large system {N = 2500) in Fig. [51 For 
optical convenience we suppress the results for the 
second and fourth sites respectively. As seen in the 
figure there may be a considerable overlap between 
the distributions corresponding to the optimal site 
and the sub-optimal ones, so that the population of 
bees may well misidentify the best site with a sig- 
nificant probability. Since the variance of the distri- 
butions shown in Fig. [H] scales as by virtue of 
the central limit theorem, this effect will generally be 
the less pronounced the larger the swarm size. For 
large swarms the individual population levels fluc- 
tuate less, resulting in sharper peaks, and a lesser 
chance of favoring a sub-optimal state. This is a spe- 
cial case of what i s known as 'Condorcet's jury the- 
orem' (see [List et al. 2 009] and references therein), 
and the effect is illustrated in Fig. where we plot 
the frequencies with which the different sites are fa- 
vored in simulations at fixed model parameters and 
system sizes. More precisely, in each simulation run, 
we check for each time step which site is the most 
populated at this moment in time, and from this 
construct the histograms shown in the figure, indi- 
cating how frequently the different sites come up as 
the most-populated one (the 'site' a = (passive 
bees) is here not taken into account). We take a sim- 
ple majority as sufficient criterion for a consensus, 
although this is admittedly a rather weak definition 
(termed 'weak consensus cr iterion' by LES) , stronger 
conditions are discussed in [List et al. 2009j | . As seen 
in the figure site a = 5 comes up with the highest 
rate for the model parameters chosen in the figure, 
but there is also a considerable probability for the 
sub-optimal sites to have the most bees dancing for 
them. From our simulations we conclude that this ef- 
fect is relevant mostly at small independence (/x close 
to unity), and tends to become stronger the smaller 
the interdependence parameter A, and/or the smaller 
the swarm. 

These finite-size stochastic effects can actually be 
understood analytically, based on an approach using 
the so-called van Kampen expansion in the inverse 
system size. In the remainder of this section we will 
describe some of the intermediate steps, but will not 
go through the derivation in all detail, as the neces- 
sary mathematics is readily available in the literature 




FIG. 5: The population levels na/N are stochastic quan- 
tities for finite systems, we here show their distributions 
for sites a = 1, a = 3 and a = 5 for fixed model pa- 
rameters /X = 0.95, A — 0.4. The system-size is A*' = 2500. 
Measurements in simulations are taken after an equilibra- 
tion period of IQQQO time-steps, and results are averaged 
over 20 independent runs. Solid lines are the theoretical 
predictions derived from Eqs. (|15|l . i.e. they are Gaussian 
distributions of mean and with variances Aaa ■ 



van Kampen 199^ . As a starting point we separate 
off stochastic fluctuations from the deterministic so- 
hition f(t) of Eqs. ©: 



(10) 



anticipating that the magnitude of deviations from 
the deterministic system scales as N~'^/'^. O ne then 
proceeds as described in [van Kampen 1993 | , and ex- 
pands the master equation ^ in powers of N~^^'^. In 
lowest order one recovers the deterministic dynamics 
of Eqs. ([7]), and in next-to-leading order one has 



(11) 



in the stationary regime, where a fixed-point of the 
mean-field trajectory has been assumed. The matrix 
J is the fc X fc Jacobian of the deterministic dynamics 
at this fixed point, i.e. one has 

Ja/3 = - [(1 - A)7r„ + A/*] + [A(l - p*) - r„] Sc^p 

(12) 

for a, /3 € {1, . . . , k}. The quantity r) in (fTTj) is Gaus- 
sian white noise of zero mean and with the following 
diagonal covariance matrix between components 



(13) 



where 




1 2 3 4 5 



1 2 3 4 5 



FIG. 6: Probabilities for the different sites a to be the 
most-populated one. Filled bars show results from sim- 
ulations for A'^ = 2500, averaged over 20 samples, open 
green bars are for A'^ = 400 (averaged over 50 samples); 
the measurement period starts after 10000 time-steps in 
the simulations. The stars and crosses show the estimates 
from the analytical theory in the first order of the system- 
size expansion, as described in the text. 



These results can also be read off from 
the general expre ssions given for example in 
Boland et al. 2008l |. Equation pl|) describes 
an multi-comp onent Ornstein-Uhlenbeck process 
[Gardiner 20"09j . and the co-variance matrix of the 
random variables £,a{t) hence has a relatively simple 
time-evolution. In partic ular if we defi ne Aa/3(t) = 



{S,a{t)£,(}{t)) , then one has Risken 1996l | 



(l-p*)[(l-A)7r„+A/:]+r„/:. (14) 



Aq/3 — ^ [Ja-yA^^ -I- J/3^A^a] -I- Fq,^, (15) 
7=1 

for a,/3 G {l,...,fc}. Denoting the fixed-point of 
these equations by A one has asymptotically 

= AAexp (y\^^W^'\p^^ (16) 

for the distribution of the population vector n. M 
is a normalization factor. Within our approximation 
(consisting in the truncation of the van-Kampen ex- 
pansion after the sub-leading order) this distribution 
is Gaussian, with a co-variance matrix A. It is here 
appropriate to point out that, within this approxi- 
mation, the expression in Eq. (fTC|) characterizes the 
statistics of the system at finite sizes fully, as both 
the variances of the individual components of n and 
their correlations are taken into account. This allows 
one to estimate mathematically how often each of the 
sites will be the 'winner' of the quorum decision, for 
example, using the weak consensus criterion of LES. 

In order to generate these semi-analytical estimates 
characterizing the finite-size effects of the system, we 
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have first integrated the master equation to asymp- 
totic times to extract the fixed-point vector f * . From 
this one computes the Jacobian J, and the above 
matrix F. Integrating Eqs. then gives A. A 

comparison of these findings against simulations is 
shown in Fig. [5l and as seen in the figure reasonably 
good agreement is achieved so that we conclude that 
our theoretical approach gives accurate estimates for 
the variances of the {na}, at least for the model pa- 
rameters used in Fig. [3 The matrix A is then in- 
verted, and subsequently results are inserted into Eq. 
(jlSp . and numerical integration of the resulting multi- 
variate Gaussian distribution then produces analyt- 
ical estimates for the probabilities for the different 
sites a G {1, . . . , fc} to be the most-populated one. 
These semi-analytical results are compared with sim- 
ulations in Fig. and the agreement is found to 
be reasonable, especially at large system sizes. We 
would here like to stress that we have not been able 
to obtain agreement of the level shown in Figs. [Hland 
ini for all values of the model parameters. In partic- 
ular if A is large, deviations can be significant, even 
for the largest system-sizes we have tried (up to 10''). 
We attribute this to remaining finite-size corrections 
(truncating the van Kampen expansion might not be 
justified in these circumstances), to equilibration ef- 
fects, or to numerical inaccuracies in integrating the 
5-dimensional Gaussian distribution (|16p . Also note 
that the Gaussian approximation becomes inaccurate 
if any of the /* are near the upper or lower limits of 
the interval [0,1]. In simulations no Ua/N can ever 
come out negative or exceed unity, these boundary 
effects are not captured by the first-order approxima- 
tion of the van Kampen expansion. As a final note in 
this section, we point out that the stro nger consensus 
criterion discussed in [List et al. 2009l| can be studied 
based on Eq. (fTH|) as well, as this equation provides 
a Gaussian estimate for the joint distribution of all k 
population levels. We have however not tried to do 
so here. 



IV. SPATIAL MODEL 
A. Model definitions and simulation results 



i.e. we introduce /^(x, t) as the fraction of bees on 
nearest neighbor lattice sites of x dancing for site a. 
More precisely we have /a(x,i) = \ Y.yr^^^a,sy(t), 
where the sum over y extends over the four nearest 
neighbors of x. A bee located at lattice site x then 
takes into account the /q;(x, f), a = 1,...,A: when 
it considers to commence a dance. The rate with 
which it starts dancing for site a is then given by 
Pc((x, i) = (1 — A)7rQ + A/a(x, i). Apart from this 
modification the dynamics proceeds exactly as in the 
non-spatial model. 

One should point out that we are here departing 
from a realistic modelling of the decision making of 
bees, as space in our model is isotropic and invari- 
ant against translations. No special point in space 
is designated as the entrance into the hive (where in- 
coming bees would arrive when they return from site- 
inspection flights). However, as discussed in the con- 
clusions, the dynamics proposed by LES can be un- 
derstood as a more general model for decision making, 
applicable also in contexts different from the nest-site 
choice of bees, so that it is worthwhile investigating 
the spatial extension discussed in the present section. 

Results from simulations are shown in Fig. [71 
When the simulations are initiated, all bees are pas- 
sive, i.e. not dancing. As seen in the figure non-trivial 
large-scale patterns emerge for sufficiently large in- 
terdependence, see especially the configurations cor- 
responding to A = 0.9 and A = 0.99 at intermediate 
times t ~ 100. Inspecting the system at later times 
(right-most snapshots in Fig. [7]) shows that these pat- 
terns are of a transient nature only, and a coarsening 
dynamics to a broad consensus on the best nest-site 
emerges in the spatial system. Interestingly, tran- 
sient structures do depend on the interdependence 
parameter A: strong interdependence promotes large 
domains of agreement on one of the different nest 
sites, whereas at small A no such large-scale areas of 
consensus are found. We should here also point out 
that these patterns do seem to depend on the initial 
conditions chosen, e.g. if we initiate the system from 
a state where all bees are active and dancing for ran- 
domly chosen sites, then the formation of domains 
seems to be far less pronounced. 



In this section we consider a spatial extension of the 
simplified model. To this end we place the dancing 
or resting bees on a square lattice of lateral extension 
L, and apply periodic boundary conditions for sim- 
plicity. The total system is composed of N = L x L 
agents. The update dynamics is as in the simplified 
model described above, the only difference is that any 
given non-dancing bee is only affected by bees in its 
direct neighborhood when it decides to start danc- 
ing for a given site, and not by the entire popula- 
tion. Specifically, the above fractions fa{t), of bees 
dancing for site a, now carry a dependence on space. 



B. Analytical approaches 



In order to make further progress towards an ana- 
lytical description of the spatial system, we here ex- 
tend the above master equation approach, and study 
the following spatial analogue of the mean-field dy- 
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FIG. 7: Pattern formation and coarsening dynamics in the spatial model at = 0. The panels show the configuration 
of the system at A = 0.1,0.9 and A = 0.99 respectively, and at different times of the simulation runs (i.e. after 
t = 10'^, 10'', 10* sweeps). Simulations are performed on a square lattice with L x L = 100^ agents, with periodic 
boundary conditions, and started from an all-passive initial condition. The different colors indicate difi'erent states of 
the individual bees, with white squares corresponding to bees who are not dancing, and filled squares to bees dancing 
for a particular site as indicated by the different colors. Lighter colors correspond to the states with low quality, and the 
black (blue on-line) areas to bees dancing for the best site. The optical impression may sometimes not reflect the actual 
proportions of bees dancing for the actual sites, accurately, e.g. in the top row we have 27%, 31% and 31% of bees 
dancing for site 5, in the bottom row these numbers are 24%, 52% and 97% at the three different times respectively. 



namics ^ of the non-spatial system. 



(1 - A)^„ 



-fA(A^„(x,t) -|-^„(x, t)) 



(17) 



for a — l,...,/c. Here (pct(x, i) is a continuous- 
valued field defined on the lattice points of the sys- 



tem, and indicates the probability of finding the agent 
at site X in state a. The (pa{'x,t) are hence coarse- 
grained order parameters, but no continuum limit 
in space is implied. Note that Eqs. (fT7|) are cou- 



pled through p{x,t) = J2a=i'^aix,t). The quan- 
tity A(pa{x,t) is the lattice Laplacian, i.e. we have 

Aipa{x,t) = (^jJ2yr^x^a{y,t)'j - (Pa(x,t). The 

term Atpa (x, t) + ipa (x, t) in the above expression 
of the right-hand-side of Eqs. pT|) is thus equal to 
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jJ2y^x'^a{y,t), and describes the total density of 
individuals in state a in the neighbourhood of x. 

Results from a numerical integration of these equa- 
tions are shown in Fig. [HI for a case in which random 
initial conditions are chosen to introduce spatial in- 
homogeneity. We find that the prediction from the 
spatial and non-spatial deterministic dynamics do not 
differ much, and that Eqs. (fT7|) capture the qualita- 
tive behavior of the dynamics, but that no quanti- 
tative agreement is achieved. The dynamics on the 
lattice appears to converge more slowly than in the 
well-mixed system, and than predicted by Eqs. (fT7| . 
These deviations point to the observation that while 
neglecting fluctuations was justified in the well-mixed 
non-spatial system, stochasticity is relevant in the 
two-dimensional system. In order to examine the spa- 
tial structures one obtains from Eqs. pT|) we show 
a number of snapshots in Fig. [HI It should be noted 
that these were taken from random initial conditions, 
with 5% of the bees active, as explained in the figure 
caption. Homogeneous initial conditions in the noise- 
less equations (|17p do not lead to any pattern forma- 
tion, as the field trivially remains homogeneous for all 
times. Fig. [9]shows domain formation similar to that 
in Fig. [71 but the extent to which the patterns of the 
agent-based model are reproduced by the determinis- 
tic dynamics of Eqs. (fTT]) is at present still undeter- 
mined, and will be left for future work. Integrating 
Eqs. (fT7)) at A = 0.9 for example (not shown) yields 
results very similar to those at A = 0.99, whereas 
in the agent-based model systematic differences are 
found between those two cases, as seen in Fig. [Tj 
To understand the relation between the agent-based 
model, and the effective dynamics in terms of con- 
tinuous fields a detailed analysis of e.g. correlation 
lengths and other statistical features of the spatio- 
temporal patterns will be required. Based on the ob- 
servations reported in Fig. [5] it may well be the case 
that a more accurate description will make it neces- 
sary to take into account fluctuations around the de- 
terministic dynamics. We have not attempted to pur- 
sue this further, but based on existing work in voter 
model s and related spati al agent-based systems (see 



where 



e.g iDornic et a l. 20051 jReichenbach al. 2008 | 



Dall'Asta fc Galla 2008| and references therein), 



one may consider the following stochastic extension 



(1 - A)^o 



iPaix,t) = 1 - p(x,i) 



+A(A¥>„(x,t) +(p„(x, t)) 

-ra(pa{x,t) + Caix,t), (18) 

where Ca (x, t) is Gaussian noise of zero mean and 
with correlations 

(Ca(x,t)C/3(x',t')> ='^a/3Baa<5(x-x')(5(i-i'), (19) 



1 - p(x, t) 



(1 - A)^„ 



-KA(A^„(x,t) -f v'a(x, t)) 

+ra(pai'X, t). 



(20) 



This covariance matrix is an analogu e of l|14 | ) and 
can be obtained following [Dall'Asta fc Galla 2008l |. 
Due to the multiplicative noise, it is difficult to inte- 
grate equations p8ll9l20p numerically, so that we do 
not present any further results here, but only men- 
tion that one approa ,ch might b e to extend recently 
developed methods |Moro 2004 iDornic et al. 20051 1 
to the case of multiple fields. We would like to 
stress that ()17|18|19|20p are not to be understood 
as a rigorously valid description of the spatial sys- 
tem. We feel that further more detailed considera- 
tions of the continuous description of the model are 
beyond the scope of the present paper, but would like 
to mention that several approaches could be consid- 
ered in future work: explicit diffusion (particle ex- 
change between neighbouring sites) could be taken 
into account, and a cont inuum limit could be de- 
rived along the lines of [Reichenbach et al. 20081. us - 
ing a Kra mers-Moyal expansion (see Gardiner 2009l | , 
and also 'Traulse n et al. 2 006'|). Alternatively sys- 
tematic expansions for spatial models with a large 
nui nber of individuals per lattice site can be found 



Lugo fc McKane 2008l | , where a van Kani pen 



system-size expansion is used [van Kampen 1992 1 



V. CONCLUSIONS 

The model of the nest-site choice dynamics in 
swarms of honeybees by LES is a beautiful example 
of a biologically inspired agent-based model, which 
can be solved with techniques originally developed 
in statistical physics. As we have shown, only mod- 
est modifications to the dynamics proposed by LES 
are required to make the model analytically tractable. 
The reduced model is Markovian in a relatively simple 
configuration space, and we have used a master equa- 
tion approach to derive a deterministic continuous- 
time description in terms of a set of coupled differen- 
tial equations. These can be integrated numerically, 
and alternatively their fixed points can be obtained 
as a the solution of a small set of quadratic equa- 
tions. The reduction may limit the degree of real- 
ism in the model, but as we have shown the simpler 
model reproduces most, if not all features of the orig- 
inal dynamics. In particular it confirms that in order 
for the population of bees to accurately identify the 
best site for a future nest, both independence in the 
assessment of the quality of the nest sites and inter- 
dependence between the bees is necessary. Based on 
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FIG. 8: Fraction of bees dancing for the best and second- 
best sites respectively (upper curves resp. lower curves). 
The smooth solid lines show results from an integration of 
the deterministic dynamics of Eqs. p7p . the dashed lines 
are from an integration of the non-spatial deterministic 
dynamics given by Eqs. (Ujl. Noisy lines show three inde- 
pendent simulation runs of the spatial agent-based model 
for A'^ = 100^ agents. Markers are from simulations of the 
non-spatial model (TV — 10*, average over 10 samples). 
Model parameters are ^ = 0, A = 0.9, in the initial state 
all bees are active and dance for randomly chosen sites. 

our analytical solutions we are able to characterize 
the stationary state of the system in the entire pa- 
rameter space spanned by the variables A and /i rep- 
resenting the degree of interdependence and the de- 
gree of independence respectively, without the need 
to perform extensive agent-based simulations. While 
the analytical description is mostly on the determin- 
istic level, a systematic approach to finite-size effects 
is feasible based on an expansion in the inverse sys- 
tem size. This allows one, within the limitations dis- 
cussed above, to obtain estimates for the frequency 
with which finite swarms of bees will reach a consen- 
sus for each of the potential nest-sites, including the 
sub-optimal ones. 

We have also extended the model to a spatial 
arrangement on a simple square lattice, and have 
demonstrated in computer simulations that the in- 
teracting agent model gives rise to non-trivial tran- 
sient patterns before a coarsening dynamics towards 
a consensus on the best nest-site sets in. These pat- 
terns are observed only at moderate to high inter- 
dependence of the bees. Finally, we have discussed 
a simple set of reaction-diffusion like equation de- 
rived from the master equation of the spatial system. 
These equations capture some of the characteristics of 



the spatial patterns formed in the agent-based model, 
but further work is required to determine whether 
they faithfully capture statistical features of these 
patterns. It is important to note that we have not 
taken into account explicit diffusion though, as no 
hopping of agents or particle exchange is considered 
in our model. A more detailed analysis of the effects 
of multiplicative noise on the field equations is a fur- 
ther interesting point, which needs a more detailed 
investigation. 

We would finally like to comment on possibilities 
to account for uncertainties in the quality assessment 
as LES do in their original model (parametrized by 
a). One way of doing this may be to formally increase 
the number of sites in our simplified model, i.e. to in- 
troduce several 'copies' of each site. E.g. for a given 
site a, one could introduce formal copies a^^\ a^^\ ... 
and then draw formal quality factors qad) , 9q(2) , . . . 
from some distribution centered around the true qual- 
ity of site a, and with a standard deviation given by 
the uncertainty parameter a. This would mathemati- 
cally increase the number of 'sites' in the model, even 
though all a'^^\ j = 1, 2, . . . would refer to the same 
real-world site, i.e. the total number of bees dancing 
for site a would be computed as n^°* = n^u) . 

Other points of further research might address the 
model in spatial arrangements different from a sim- 
ple two-dimensional regular lattice. One may here 
for example study the effects of the topology of the 
network of interactions on the convergence and con- 
sensus properties. In this context the model proposed 
by LES and the simplified dynamics discussed in the 
present paper could probably be interpreted more 
generally as a model for the interplay of indepen- 
dence and interdependence in the decision making of 
groups, where the specific application may not neces- 
sarily be limited to the nest-site choice of honeybees. 
LES give an example towards the end of their pa- 
per, relating to the choice of restaurants by humans, 
and similarly many other situations, where individu- 
als have to make strategic choices, for example in the 
social sciences and in game theory may potentially be 
modeled using similar dynamics. The analytical ap- 
proaches developed in the present paper may here be 
useful for to characterize the outcome of the decision 
making processes in such models. 
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